Intersections of Adversity & Neurodiversity: Adverse Childhood Experiences’ Association to Mental Health & the Buffering Role of Flourishing
This study extends analyses conducted by Bethell et al. (2023), who explored the role of family resilience and connection in promoting flourishing among U.S. children facing adversity. Our study applies a similar framework but focuses on neurodivergent children, examining how adverse childhood experiences (ACEs) influence mental health outcomes and how flourishing behaviors—such as curiosity, emotional control, and task persistence—may buffer these effects. This study was conducted by the Slopen Laboratory at the Harvard T.H. Chan School of Public Health.
Author: Adrián Medina
Date: November 24, 2024
Project Overview
This repository contains files relevant to a study investigating the relationship between adverse childhood experiences (ACEs) and mental health outcomes in children with neurodevelopmental disorders (NDDs), such as Autism Spectrum Disorder (ASD) and ADHD. The study specifically examines the moderating role of child flourishing behaviors—curiosity in learning, emotional control, and task persistence—on the effects of ACEs on mental health outcomes.
Study Design:
The study utilizes data from the National Survey of Children’s Health (NSCH), focusing on children aged 6-17 with reported neurodevelopmental issues (N = 44,776, MAge = 12.2). This cross-sectional analysis examines ACE exposure and its impact on three primary mental health outcomes: anxiety, depression, and behavioral issues. The moderating role of child flourishing behaviors is assessed using interaction terms in logistic regression models.
Data Collection:
Key variables collected include:
Adverse Childhood Experiences (ACEs): Derived from parent-reported experiences such as exposure to violence, household dysfunction, and discrimination.
Mental Health Outcomes: Parent-reported diagnoses of anxiety, depression, and behavioral problems, confirmed by health care providers or educators.
Child Flourishing Metrics: Based on parent-reported behaviors such as curiosity in learning, emotional control, and task persistence, which were used to develop a Child Flourishing Index (CFI).
Analytic Approach:
The primary analytic approach involved logistic regression models to assess the association between ACEs and mental health outcomes, with interaction terms to explore the moderating effect of child flourishing. Key covariates include age, sex, race/ethnicity, and socioeconomic status.
Goals:
The study aims to:
Examine the dose-response relationship between the number of ACEs and risks of mental health challenges in neurodivergent children.
Identify the protective role of child flourishing behaviors, offering insights into resilience-promoting interventions for this vulnerable population.
Tip
For detailed information on the dataset variables, refer to the NSCH Data Dictionary.
Set up the R environment by configuring the CRAN repository, installing essential packages, setting the base path, and loading the primary dataset into the local environment.
Expand Code
# Set CRAN repository for consistent package installationoptions(repos =c(CRAN ="http://cran.rstudio.com/"))# Install and load the pacman package for efficient package managementif (!require(pacman)) install.packages("pacman")library(pacman)# Use p_load function to install (if necessary) and load packagesp_load(dplyr, tidyr, tidyverse, knitr, ggplot2, Hmisc, broom, stats, kableExtra, webshot2, sjPlot, margins, gtsummary, ggstats, broom.helpers, patchwork, sjlabelled, sjmisc, vtable, DiagrammeR)# Specify the 'base_path' where you can find your data files, ASSUMING they're in the same directory, & set it as WDbase_path <-"~/GitHub/Adversity-Neurodiversity/data_files"setwd(base_path)# Load primary data fileneurodivergent_data <-read_csv("neurodivergent_data.csv")
1
See details under the Data Extraction & Filtering (Archived Code) callout below!
Data Extraction & Filtering (Archived Code)
The file being used for these analyses is a subset of a “master” file, Data2_2016to2022.csv, containing the compiled contents of the data releases from the National Survey on Child Health’s 2016-2022 cycles. The master file is nearly 500 MB (0.5 GB) in size, making it cumbersome to process in real-time. For efficiency and ease of analysis, I am using a smaller subset that contains only the relevant data. This subset was generated using the data extraction and filtering process detailed in the archived code below, in an effort to maintain transparency and ensure reproducibility of the workflow.
Expand Code
# Define the path to the datasetdata_path <-"~/Downloads/Data2_2016to2022.csv"Data2_2016to2022 <-read.csv(data_path)# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent statusneurodivergent_data <- Data2_2016to2022 %>%select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis, k2q35a, k2q35c, k2q38a, k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, k2q37a, k2q37c, downsyn, downsyn_desc, k2q31a, k2q31c, k2q61a, cerpals_desc, k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, k7q85_r, k7q84_r) %>%# Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.# Using a list of diagnoses provided by NSCH to define neurodivergent individuals.mutate(neurodivergent =if_else(k2q35a ==1| k2q38a ==1| k2q36a ==1| k2q60a ==1| k2q37a ==1| k2q30a ==1| downsyn ==1| k2q31a ==1| k2q61a ==1, 1, 0)) %>%# Filter data to include only individuals aged 6 or above and identified as neurodivergentfilter(sc_age_years >=6, neurodivergent ==1)# Resulting filtered data to be used for further analysiswrite.csv(neurodivergent_data, "neurodivergent_data.csv")
Analytic Data Preparation
Recode analytic variables and covariates, including ACE counts, Childhood Flourishing Index (CFI), binary outcomes (anxiety, depression, behavioral issues), education, income, sex, and race/ethnicity variables into categorical factors for analysis.
Expand Code
# Recode 'ace1' to dichotomous variable based on analytic specificationsneurodivergent_data$ace1_recode <-ifelse(neurodivergent_data$ace1 %in%c(2, 3, 4), 1, ifelse(neurodivergent_data$ace1 ==1, 0, NA))# Recode ACE counts and calculate Childhood Flourishing Index (CFI)neurodivergent_data <- neurodivergent_data %>%mutate(# Total ACE count for each child excluding missing valuesace_total =rowSums(select(., ace1_recode, ace3:ace12) ==1, na.rm =TRUE),# Categorize total ACEs for further analysisace_counts =cut(ace_total, breaks =c(-1, 0, 1, 3, Inf), labels =c('0 ACEs', '1 ACE', '2-3 ACEs', '4+ ACEs'), right =TRUE),ace_counts =factor(ace_counts, levels =c("0 ACEs", "1 ACE", "2-3 ACEs", "4+ ACEs")), # Recreate 'ace_counts' factor# Calculate the Childhood Flourishing Index (CFI)CFI =rowSums(select(., k6q71_r, k7q85_r, k7q84_r) ==1, na.rm =TRUE),# Recode CFI into dichotomous variableCFI_dich =case_when( CFI ==3~"Flourishing", CFI %in%0:2~"Not Flourishing",TRUE~NA_character_ ),CFI_dich =factor(CFI_dich, levels =c("Not Flourishing", "Flourishing")) # Recreate 'CFI_dich' factor )# Recode binary outcomes for anxiety, depression, and behavioral issues into factorsneurodivergent_data <- neurodivergent_data %>%mutate(# Recode Anxiety (k2q33b)Anxiety =factor(k2q33b, levels =c(2, 1), labels =c("No", "Yes")),# Recode Depression (k2q32b)Depression =factor(k2q32b, levels =c(2, 1), labels =c("No", "Yes")),# Recode Behavioral Problems (k2q34b)Behavioral_Problems =factor(k2q34b, levels =c(2, 1), labels =c("No", "Yes")) )# Recode education and income categoriesneurodivergent_data <- neurodivergent_data %>%mutate(# Recode education into categorical factorshighgrade_tvis_cat =case_when( higrade_tvis ==1~"Less Than High School", higrade_tvis ==2~"High School/Vocational", higrade_tvis ==3~"Some College/Associate Degree", higrade_tvis ==4~"College Degree/Higher",TRUE~NA_character_ ),highgrade_tvis_cat =factor(highgrade_tvis_cat, levels =c("Less Than High School", "High School/Vocational", "Some College/Associate Degree", "College Degree/Higher")),# Recode federal poverty level categories into factorsfpl_cat =case_when( fpl %in%c(50:99) | fpl_mean <100~"<100%", fpl %in%c(100:199) | fpl_mean <200~"100%-199%", fpl %in%c(200:299) | fpl_mean <300~"200%-299%", fpl %in%c(300:399) | fpl_mean <400~"300%-399%", fpl %in%c(400:999) | fpl_mean <999~"≥400%",TRUE~NA_character_ ),fpl_cat =factor(fpl_cat, levels =c("<100%", "100%-199%", "200%-299%", "300%-399%", "≥400%")) )# Recode sex and race categoriesneurodivergent_data <- neurodivergent_data %>%mutate(# Recode sex into categorical factorssc_sex_cat =case_when( sc_sex ==1~"Male", sc_sex ==2~"Female",TRUE~NA_character_ ),sc_sex_cat =factor(sc_sex_cat, levels =c("Male", "Female")),# Recode race and ethnicityrace =case_when( sc_race_r ==1~"White", sc_race_r ==2~"Black", sc_race_r ==3~"American Indian or Alaska Native", sc_race_r %in%4:5~"Asian, Native Hawaiian, or Pacific Islander", sc_race_r %in%6:7~"Other or Mixed Race",TRUE~NA_character_ ),ethnicity =case_when( sc_hispanic_r ==1~"Hispanic/Latino", sc_hispanic_r ==2~"Non-Hispanic/Latino",TRUE~NA_character_ ),# Recreate 'sc_race_cat' from race and ethnicity combinationssc_race_cat =case_when( race %in%c("White") & ethnicity =="Non-Hispanic/Latino"~"Non-Hispanic White", race %in%c("Black") & ethnicity =="Non-Hispanic/Latino"~"Non-Hispanic Black or African American", race %in%c("American Indian or Alaska Native") & ethnicity =="Non-Hispanic/Latino"~"Non-Hispanic American Indian or Alaska Native", race %in%c("Asian, Native Hawaiian, or Pacific Islander") & ethnicity =="Non-Hispanic/Latino"~"Non-Hispanic Asian, Native Hawaiian, or Pacific Islander", race %in%c("Other or Mixed Race") & ethnicity =="Non-Hispanic/Latino"~"Non-Hispanic Other or Mixed Race", ethnicity =="Hispanic/Latino"~"Hispanic/Latino of Any Race",TRUE~NA_character_ ),sc_race_cat =factor(sc_race_cat, levels =c("Non-Hispanic White", "Non-Hispanic Black or African American", "Non-Hispanic American Indian or Alaska Native", "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander", "Non-Hispanic Other or Mixed Race", "Hispanic/Latino of Any Race")) )
1
ace1 is originally coded by NSCH as a frequency measure, but it’s needed as indicative variable of presence like the other ace# variables.
2
Part of this includes the calculating of the Childhood Flourishing Index (CFI) based on Bethell et al. (2019) criteria
Frequencies & Descriptive Statistics
Calculate the prevalence of neurodevelopmental disorder (NDD) diagnoses by defining the relevant variables, computing the number of total respondents for each variable, and determining the percentage of diagnoses. Display the results in a table with a table note explaining the presentation format.
Expand Code
# Calculate prevalence (%) of NDD diagnoses among subject samplevariables <-c("k2q31a", "k2q35a", "k2q38a", "k2q30a", "k2q36a", "k2q60a", "k2q37a", "downsyn", "k2q61a")variable_names <-c("ADHD", "Autism/ASD", "Tourette’s Syndrome", "Learning Disability", "Development Delay", "Intellectual Disability", "Speech/Other Language Disorder", "Down Syndrome", "Cerebral Palsy")# Pre-calculate total respondents to avoid repeated computationtotal_respondents <-colSums(!is.na(neurodivergent_data[variables]))# Vectorized prevalence calculationprevalences <-colSums(neurodivergent_data[variables] ==1, na.rm =TRUE) / total_respondents *100prevalence_table <-data.frame(Diagnosis = variable_names, Prevalence = prevalences)# Display prevalence data in a tableprevalence_table %>% knitr::kable("html", caption ="Prevalence of Neurodivergent Conditions") %>% kableExtra::footnote(general ="Prevalence values are presented as percentages for neurodevelopmental disorder diagnosis.")
Prevalence of Neurodivergent Conditions
Diagnosis
Prevalence
k2q31a
ADHD
58.9175975
k2q35a
Autism/ASD
15.2625563
k2q38a
Tourette’s Syndrome
1.6302642
k2q30a
Learning Disability
39.6132560
k2q36a
Development Delay
31.7615419
k2q60a
Intellectual Disability
5.6693867
k2q37a
Speech/Other Language Disorder
36.2469766
downsyn
Down Syndrome
0.9816886
k2q61a
Cerebral Palsy
1.4784403
Note:
Prevalence values are presented as percentages for neurodevelopmental disorder diagnosis.
Summarize key demographic characteristics of the neurodivergent sample, including descriptive statistics for age, educational attainment, sex, socioeconomic status, and race. Display the results in a table and visualize distributions using raincloud and bar plots for each demographic variable.
Expand Code
# Calculate descriptives of age for neurodivergent participantsneurodiv_age_summary <- neurodivergent_data %>% dplyr::summarise(score_mean =mean(sc_age_years, na.rm =TRUE), # Meann =n(), # Sample sizesem =sd(sc_age_years, na.rm =TRUE) /sqrt(n()), # Standard error of the meansd =sd(sc_age_years, na.rm =TRUE), # Standard deviationmin_age =min(sc_age_years, na.rm =TRUE), # Minimumq1_age =quantile(sc_age_years, 0.25, na.rm =TRUE), # 1st Quartilemedian_age =median(sc_age_years, na.rm =TRUE), # Medianq3_age =quantile(sc_age_years, 0.75, na.rm =TRUE), # 3rd Quartilemax_age =max(sc_age_years, na.rm =TRUE), # MaximumIQR_age =quantile(sc_age_years, 0.75, na.rm =TRUE) -quantile(sc_age_years, 0.25, na.rm =TRUE), # IQR (Q3 - Q1).groups ='drop' )# Create a neat table using kableneurodiv_age_summary %>% knitr::kable("html", col.names =c("Mean", "Sample Size", "SEM", "SD", "Min", "Q1", "Median", "Q3", "Max", "IQR"), caption ="Descriptive Statistics of Age for Neurodivergent Participants")
Descriptive Statistics of Age for Neurodivergent Participants
Mean
Sample Size
SEM
SD
Min
Q1
Median
Q3
Max
IQR
12.17523
44776
0.0159532
3.375748
6
9
12
15
17
6
Expand Code
# Raincloud plot of ages for neurodivergent participants with colorggplot(neurodivergent_data, aes(x =1, y = sc_age_years)) +# Set x to 1, as there's only one group PupillometryR::geom_flat_violin(adjust =1.5, trim =FALSE, alpha =0.5, color =NA, fill ="deeppink", position =position_nudge(x =0.1, y =0)) +geom_point(alpha =0.1, position =position_jitter(width = .05), size = .10, shape =20, color ="deeppink") +geom_boxplot(outlier.shape =NA, alpha =0, width =0.1, position =position_dodge(width =0.3), colour ="black") +geom_point(data = neurodiv_age_summary, aes(x =1, y = score_mean), shape =18, color ="deeppink", size =3, position =position_nudge(x =0.1)) +geom_errorbar(data = neurodiv_age_summary, aes(x =1, y = score_mean, ymin = score_mean - sem, ymax = score_mean + sem), width =0.05, color ="deeppink", position =position_nudge(x =0.1)) +labs(title ="Age Distribution of Neurodivergent Participants", y ="Age (Years)", x ="") +# No x-axis label since there's only one groupscale_y_continuous(breaks =seq(5, 18, by =1), limits =c(5, 18)) +# Set y-axis range and incrementstheme_minimal() +theme(axis.text.x =element_blank()) # Remove x-axis text
Calculate and display the number of subjects classified as ‘Flourishing’ based on a Childhood Flourishing Index (CFI) score of 3. Visualize the distribution of children by CFI score using bar plots for both the raw CFI scores and the dichotomized CFI variable.
Expand Code
# Count the number of subjects with a CFI score of 3 using a logical condition, as this is considered 'Flourishing.'num_subjects_cfi_3 <-sum(neurodivergent_data$CFI ==3, na.rm =TRUE)cat("Number of Subjects Flourishing (CFI Score of 3):", num_subjects_cfi_3, "\n")
Number of Subjects Flourishing (CFI Score of 3): 3265
Expand Code
# Bar plot of CFI counts categoriesggplot(neurodivergent_data, aes(x =factor(CFI), fill =factor(CFI))) +geom_bar(show.legend =FALSE) +geom_text(stat ='count', aes(label = ..count..), position =position_stack(vjust =0.5), color ="white") +labs(x ="Childhood Flourishing Index (CFI) Score",y ="Number of Children",fill ="CFI Score",title ="Distribution of Children by CFI Score") +theme_minimal() +theme(axis.text.x =element_text(angle =0))
Expand Code
# Bar plot of CFI_dich counts categoriesggplot(neurodivergent_data, aes(x = CFI_dich)) +geom_bar(aes(fill = CFI_dich), show.legend =FALSE) +geom_text(stat ='count', aes(label = ..count..), position =position_stack(vjust =0.5), color ="white") +labs(title ="Distribution of CFI Dichotomous Variable",x ="Child Flourishing Index (Dichotomous)",y ="Count",fill ="CFI") +theme_minimal() +theme(axis.text.x =element_text(angle =0))
Construct and display a frequency table that shows the number of participants with Anxiety, Depression, and Behavioral Issues based on responses indicating “Yes” for each condition. Present the results in an HTML table format using kable.
Expand Code
# Constructing the frequency table for Anxiety, Depression, and Behavioral Issues where the response is "Yes"neurodiv_freq_table <- neurodivergent_data %>% dplyr::summarise(Anxiety_Yes =sum(Anxiety =="Yes", na.rm =TRUE),Depression_Yes =sum(Depression =="Yes", na.rm =TRUE),Behavioral_Problems_Yes =sum(Behavioral_Problems =="Yes", na.rm =TRUE) )# Display table in a simple HTML format using kableneurodiv_freq_table %>% knitr::kable("html", col.names =c("Anxiety", "Depression", "Behavioral Problems"), caption ="Frequency of Mental Health Conditions in Neurodivergent Sample")
Frequency of Mental Health Conditions in Neurodivergent Sample
Anxiety
Depression
Behavioral Problems
13725
6289
13535
Visualize the distribution of ACE counts among children using a bar plot and calculate cross-tabulations for the number and percentage of participants with Anxiety, Depression, Behavioral Issues, and Flourishing status across ACE categories. Format the table to display values as “n (percentage)” for each ACE category and include a footnote explaining the presentation format.
Expand Code
# Bar plot of ACE counts categoriesggplot(neurodivergent_data, aes(x =factor(ace_counts), fill =factor(ace_counts))) +geom_bar(show.legend =FALSE) +geom_text(stat ='count', aes(label = ..count..), position =position_stack(vjust =0.5), color ="white") +labs(x ="ACE Counts", y ="Number of Children", fill ="ACE Counts",title ="Distribution of Children by ACE Counts") +theme_minimal() +theme(axis.text.x =element_text(angle =0))
1
Cross-tabulations (or contingency tables) summarize the relationship between two or more categorical variables. In this case, I’m exploring how the mental health conditions (Anxiety, Depression, Behavioral Issues) and Child Flourishing Index (CFI) are distributed across different levels of ACEs.
Expand Code
# Define the total number of respondents for each ACE categoryace_totals <-c("0 ACEs"=12431, # Total respondents with 0 ACEs"1 ACE"=13433, # Total respondents with 1 ACE"2-3 ACEs"=12180, # Total respondents with 2-3 ACEs"4+ ACEs"=6732# Total respondents with 4+ ACEs)# Compute the cross-tabulations, filtering out NA values and excluding "No" responses for Anxiety, Depression, and Behavioral Issuescross_tabs <- neurodivergent_data %>%select(Anxiety, Depression, Behavioral_Problems, CFI_dich, ace_counts) %>%pivot_longer(cols =c(Anxiety, Depression, Behavioral_Problems, CFI_dich), names_to ="Condition", values_to ="Status") %>%filter(!is.na(Status), !is.na(ace_counts), !(Condition %in%c("Anxiety", "Depression", "Behavioral_Problems") & Status =="No")) %>%# Filter out NA values and "No" responsesgroup_by(Condition, Status, ace_counts) %>%summarise(n =n(), .groups ="drop") %>%group_by(Condition, Status) %>%# Calculate percentage by dividing n by the total number of respondents for each ACE category and create a string in the format n (percent)summarise(`0 ACEs`=if_else(any(ace_counts =="0 ACEs"), paste0(n[ace_counts =="0 ACEs"], " (", round(n[ace_counts =="0 ACEs"] / ace_totals["0 ACEs"] *100, 1), ")"), NA_character_),`1 ACE`=if_else(any(ace_counts =="1 ACE"), paste0(n[ace_counts =="1 ACE"], " (", round(n[ace_counts =="1 ACE"] / ace_totals["1 ACE"] *100, 1), ")"), NA_character_),`2-3 ACEs`=if_else(any(ace_counts =="2-3 ACEs"), paste0(n[ace_counts =="2-3 ACEs"], " (", round(n[ace_counts =="2-3 ACEs"] / ace_totals["2-3 ACEs"] *100, 1), ")"), NA_character_),`4+ ACEs`=if_else(any(ace_counts =="4+ ACEs"), paste0(n[ace_counts =="4+ ACEs"], " (", round(n[ace_counts =="4+ ACEs"] / ace_totals["4+ ACEs"] *100, 1), ")"), NA_character_) )# Format the table for presentation using 'kable' with the footnote spanning the full widthcross_tabs %>% knitr::kable("html", caption ="Frequency and Percentage of Psychopathology and Flourishing across ACEs Categories") %>% kableExtra::kable_styling(full_width =TRUE) %>%# Ensures the table spans the full width kableExtra::footnote(general ="Values are presented as n (percentage) for each ACEs category. Cross-tabulations show the number and percentage of individuals with a specific mental health condition or flourishing status across different ACE categories.",footnote_as_chunk =TRUE) # Formats the footnote to span full width
Frequency and Percentage of Psychopathology and Flourishing across ACEs Categories
Condition
Status
0 ACEs
1 ACE
2-3 ACEs
4+ ACEs
Anxiety
Yes
2695 (21.7)
3549 (26.4)
4257 (35)
3224 (47.9)
Behavioral_Problems
Yes
2530 (20.4)
3561 (26.5)
4189 (34.4)
3255 (48.4)
CFI_dich
Not Flourishing
11140 (89.6)
12358 (92)
11492 (94.4)
6521 (96.9)
CFI_dich
Flourishing
1291 (10.4)
1075 (8)
688 (5.6)
211 (3.1)
Depression
Yes
804 (6.5)
1222 (9.1)
2112 (17.3)
2151 (32)
Note: Values are presented as n (percentage) for each ACEs category. Cross-tabulations show the number and percentage of individuals with a specific mental health condition or flourishing status across different ACE categories.
Inferential Statistics
Run three separate logistic regression models to examine the effects of ACE counts, race, sex, socioeconomic status, caregiver education, and child age on Anxiety, Depression, and Behavioral Issues outcomes. Present the exponentiated coefficients (odds ratios) and bold significant p-values (p < 0.05). Merge the summary tables of all three models into one, with each model labeled appropriately for clear comparison.
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander
0.63
0.43, 0.95
0.021
1.34
0.83, 2.30
0.3
1.14
0.84, 1.58
0.4
Non-Hispanic Other or Mixed Race
0.89
0.71, 1.14
0.4
1.01
0.80, 1.29
>0.9
0.93
0.78, 1.10
0.4
Hispanic/Latino of Any Race
0.71
0.59, 0.85
<0.001
0.84
0.69, 1.02
0.070
0.90
0.78, 1.05
0.2
Child Sex
Male
—
—
—
—
—
—
Female
1.34
1.18, 1.51
<0.001
1.42
1.25, 1.62
<0.001
1.21
1.09, 1.35
<0.001
Federal Poverty Level
<100%
—
—
—
—
—
—
100%-199%
0.94
0.75, 1.17
0.6
0.89
0.71, 1.12
0.3
0.91
0.77, 1.08
0.3
200%-299%
1.03
0.82, 1.31
0.8
0.84
0.66, 1.06
0.14
0.83
0.70, 0.99
0.043
300%-399%
0.82
0.64, 1.04
0.10
0.79
0.61, 1.01
0.062
0.68
0.57, 0.82
<0.001
≥400%
0.91
0.72, 1.14
0.4
0.69
0.54, 0.87
0.002
0.71
0.59, 0.84
<0.001
Caregiver Education
Less Than High School
—
—
—
—
—
—
High School/Vocational
1.18
0.77, 1.77
0.4
1.10
0.71, 1.64
0.7
1.01
0.73, 1.37
>0.9
Some College/Associate Degree
1.18
0.77, 1.74
0.4
1.13
0.74, 1.68
0.6
0.99
0.72, 1.34
>0.9
College Degree/Higher
1.30
0.85, 1.93
0.2
1.03
0.67, 1.53
>0.9
0.87
0.63, 1.17
0.4
Child Age
0.97
0.95, 0.99
<0.001
0.96
0.93, 0.98
<0.001
0.86
0.84, 0.87
<0.001
1 OR = Odds Ratio, CI = Confidence Interval
Generate a combined plot showing the odds ratios from three separate logistic regression models (Anxiety, Depression, and Behavioral Issues). Configure the figure dimensions to 8x8, with confidence intervals (95%) and p-values displayed. Use different markers and colors for each model, set axis limits between 0.2 and 5, and include a legend labeling each model accordingly.
Expand Code
# Combine the plots into a list vectorall.models <-list()all.models[[1]] <- model_anxall.models[[2]] <- model_depall.models[[3]] <- model_beh# Plot the combined versionplot_models(all.models, axis.lim =c(0.2, 5), dot.size =1, show.p =TRUE, ci.lvl =0.95, title ="Odds Ratios by Model", legend.title ="Model", m.labels =c("Anxiety", "Depression", "Behavioral Problems"))
Run logistic regression models that incorporate interaction terms between the ACE counts and the Child Flourishing Index (CFI_dich) for three different mental health outcomes (Anxiety, Depression, Behavioral Issues). Then, create individual summary tables for each model (without exponentiating the coefficients), highlighting statistically significant p-values. Finally, merge the tables into a single output, displaying all models side by side with appropriate labels and formatting.
Expand Code
# Logistic regression model with CFI_dich as a moderator for Anxietymodel_anx_mod <-glm(Anxiety ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)# Create a summary table for the Anixety modeltbl_anx_mod <-tbl_regression(model_anx_mod, exponentiate =FALSE, label =list( ace_counts ~"ACEs", CFI_dich ~"Child Flourishing Index", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Logistic regression model with CFI_dich as a moderator for Depressionmodel_dep_mod <-glm(Depression ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)# Create a summary table for the Depression modeltbl_dep_mod <-tbl_regression(model_dep_mod, exponentiate =FALSE, label =list( ace_counts ~"ACEs", CFI_dich ~"Child Flourishing Index", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Logistic regression model with CFI_dich as a moderator for Behavioral Problemsmodel_beh_mod <-glm(Behavioral_Problems ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)# Create a summary table for the Behavioral Problems modeltbl_beh_mod <-tbl_regression(model_beh_mod, exponentiate =FALSE, label =list( ace_counts ~"ACEs", CFI_dich ~"Child Flourishing Index", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Merge all tables into onemod_tbl_merged <-tbl_merge(tbls =list(tbl_anx_mod, tbl_dep_mod, tbl_beh_mod),tab_spanner =c("Anixety Model", "Depression Model", "Behavioral Problems Model"))# Display the combined tablemod_tbl_merged
1
The interaction anxiety model did not yield significant findings, so we do not proceed with post-hoc analyses.
2
The interaction depression model did yield significant findinds, so we do proceed with post-hoc analyses.
3
The interaction behavioral problems model did yield significant findinds, so we do proceed with post-hoc analyses.
Characteristic
Anixety Model
Depression Model
Behavioral Problems Model
log(OR)1
95% CI1
p-value
log(OR)1
95% CI1
p-value
log(OR)1
95% CI1
p-value
ACEs
0 ACEs
—
—
—
—
—
—
1 ACE
0.26
0.09, 0.43
0.003
0.19
-0.01, 0.39
0.065
0.14
0.00, 0.28
0.043
2-3 ACEs
0.38
0.20, 0.55
<0.001
0.32
0.13, 0.51
<0.001
0.11
-0.03, 0.25
0.12
4+ ACEs
0.58
0.38, 0.78
<0.001
0.54
0.33, 0.74
<0.001
0.31
0.16, 0.46
<0.001
Child Flourishing Index
Not Flourishing
—
—
—
—
—
—
Flourishing
-0.99
-1.4, -0.50
<0.001
-1.4
-2.2, -0.61
<0.001
-0.71
-1.3, -0.13
0.013
Child Race
Non-Hispanic White
—
—
—
—
—
—
Non-Hispanic Black or African American
-0.27
-0.54, 0.01
0.050
-0.06
-0.33, 0.23
0.7
0.03
-0.15, 0.22
0.7
Non-Hispanic American Indian or Alaska Native
-0.66
-1.2, -0.07
0.018
0.28
-0.40, 1.1
0.5
-0.11
-0.57, 0.41
0.7
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander
-0.48
-0.85, -0.07
0.016
0.29
-0.19, 0.83
0.3
0.11
-0.20, 0.44
0.5
Non-Hispanic Other or Mixed Race
-0.13
-0.36, 0.12
0.3
0.01
-0.23, 0.25
>0.9
-0.07
-0.24, 0.10
0.4
Hispanic/Latino of Any Race
-0.34
-0.52, -0.16
<0.001
-0.18
-0.37, 0.02
0.075
-0.10
-0.25, 0.05
0.2
Child Sex
Male
—
—
—
—
—
—
Female
0.30
0.18, 0.42
<0.001
0.36
0.23, 0.49
<0.001
0.19
0.09, 0.30
<0.001
Federal Poverty Level
<100%
—
—
—
—
—
—
100%-199%
-0.06
-0.28, 0.16
0.6
-0.11
-0.34, 0.12
0.4
-0.11
-0.29, 0.06
0.2
200%-299%
0.03
-0.20, 0.27
0.8
-0.18
-0.42, 0.06
0.14
-0.20
-0.38, -0.02
0.026
300%-399%
-0.19
-0.43, 0.05
0.12
-0.24
-0.50, 0.01
0.062
-0.40
-0.59, -0.21
<0.001
≥400%
-0.09
-0.32, 0.14
0.5
-0.37
-0.61, -0.14
0.002
-0.37
-0.55, -0.19
<0.001
Caregiver Education
Less Than High School
—
—
—
—
—
—
High School/Vocational
0.14
-0.29, 0.55
0.5
0.09
-0.34, 0.50
0.7
-0.01
-0.34, 0.30
>0.9
Some College/Associate Degree
0.13
-0.30, 0.52
0.5
0.12
-0.30, 0.51
0.6
-0.04
-0.36, 0.27
0.8
College Degree/Higher
0.22
-0.21, 0.61
0.3
0.02
-0.41, 0.42
>0.9
-0.17
-0.49, 0.14
0.3
Child Age
-0.03
-0.05, -0.01
0.002
-0.04
-0.07, -0.02
<0.001
-0.15
-0.17, -0.14
<0.001
ACEs * Child Flourishing Index
1 ACE * Flourishing
-0.23
-0.91, 0.46
0.5
-0.13
-1.3, 1.0
0.8
-0.52
-1.3, 0.24
0.2
2-3 ACEs * Flourishing
0.10
-0.58, 0.80
0.8
1.3
0.27, 2.5
0.017
-0.54
-1.3, 0.20
0.2
4+ ACEs * Flourishing
0.54
-0.44, 1.7
0.3
1.0
-0.08, 2.2
0.077
-1.1
-2.1, -0.19
0.019
1 OR = Odds Ratio, CI = Confidence Interval
Subset the data into flourishing and non-flourishing groups. Then, run logistic regression models for Depression and Behavioral Problems separately for each group, adjusting for covariates (ACE counts, race, sex, SES, education, and age). Create summary tables for each model, with p-values highlighted for significance. Finally, merge and display the tables for both flourishing and non-flourishing samples, stratified by Depression and Behavioral Problems.
Expand Code
# Subset the flourishing sampleflourishing_sample <- neurodivergent_data %>%filter(CFI_dich =="Flourishing")# Subset the not flourishing samplenot_flourishing_sample <- neurodivergent_data %>%filter(CFI_dich =="Not Flourishing")# Stratified logistic regression modelsflourish_model_dep <-glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = flourishing_sample, family = binomial)non_flourish_model_dep <-glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = not_flourishing_sample, family = binomial)flourish_model_beh <-glm(Behavioral_Problems ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = flourishing_sample, family = binomial)non_flourish_model_beh <-glm(Behavioral_Problems ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = not_flourishing_sample, family = binomial)# Create a summary table for the Flourishing Depression modeltbl_dep_flourish <-tbl_regression(flourish_model_dep, exponentiate =TRUE, label =list( ace_counts ~"ACEs", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Create a summary table for the Non-Flourishing Depression modeltbl_dep_non_flourish <-tbl_regression(non_flourish_model_dep, exponentiate =TRUE, label =list( ace_counts ~"ACEs", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Create a summary table for the Flourishing Behavioral Problems modeltbl_beh_flourish <-tbl_regression(flourish_model_beh, exponentiate =TRUE, label =list( ace_counts ~"ACEs", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Create a summary table for the Non-Flourishing Behavioral Problems modeltbl_beh_non_flourish <-tbl_regression(non_flourish_model_beh, exponentiate =TRUE, label =list( ace_counts ~"ACEs", fpl_cat ~"Federal Poverty Level", highgrade_tvis_cat ~"Caregiver Education", sc_age_years ~"Child Age", sc_race_cat ~"Child Race", sc_sex_cat ~"Child Sex" )) |>bold_p(t =0.05)# Merge the flourishing model tables into oneflourish_tbl_combined <-tbl_merge(tbls =list(tbl_dep_flourish, tbl_beh_flourish),tab_spanner =c("Depression Model", "Behavioral Problems Model"))# Display the combined flourishing model tablesflourish_tbl_combined
Characteristic
Depression Model
Behavioral Problems Model
OR1
95% CI1
p-value
OR1
95% CI1
p-value
ACEs
0 ACEs
—
—
—
—
1 ACE
1.18
0.34, 4.11
0.8
0.71
0.32, 1.57
0.4
2-3 ACEs
7.81
2.37, 28.8
0.001
0.64
0.29, 1.38
0.3
4+ ACEs
8.99
2.36, 38.1
0.002
0.36
0.12, 1.05
0.064
Child Race
Non-Hispanic White
—
—
—
—
Non-Hispanic Black or African American
0.42
0.09, 2.09
0.3
0.80
0.29, 2.28
0.7
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander
20,679,597
0.00,
>0.9
3,087,910
0.00,
>0.9
Non-Hispanic Other or Mixed Race
0.68
0.10, 6.07
0.7
1.42
0.47, 4.94
0.6
Hispanic/Latino of Any Race
1.03
0.29, 4.22
>0.9
0.40
0.16, 0.96
0.042
Non-Hispanic American Indian or Alaska Native
0.20
0.01, 5.72
0.3
Child Sex
Male
—
—
—
—
Female
1.12
0.48, 2.63
0.8
1.72
0.88, 3.48
0.12
Federal Poverty Level
<100%
—
—
—
—
100%-199%
1.65
0.41, 7.05
0.5
0.73
0.27, 1.93
0.5
200%-299%
1.60
0.33, 8.40
0.6
0.99
0.35, 2.75
>0.9
300%-399%
0.84
0.18, 4.05
0.8
0.93
0.30, 2.87
>0.9
≥400%
2.62
0.58, 12.5
0.2
0.68
0.23, 1.92
0.5
Caregiver Education
Less Than High School
—
—
—
—
High School/Vocational
0.94
0.04, 8.79
>0.9
0.90
0.20, 3.61
0.9
Some College/Associate Degree
0.44
0.02, 3.64
0.5
1.00
0.22, 3.95
>0.9
College Degree/Higher
0.33
0.01, 2.97
0.4
0.62
0.13, 2.53
0.5
Child Age
0.91
0.75, 1.10
0.3
0.86
0.78, 0.95
0.003
1 OR = Odds Ratio, CI = Confidence Interval
Expand Code
# Merge the non-flourishing model tables into onenon_flourish_tbl_combined <-tbl_merge(tbls =list(tbl_dep_non_flourish, tbl_beh_non_flourish),tab_spanner =c("Depression Model", "Behavioral Problems Model"))# Display the combined non-flourishing model tablesnon_flourish_tbl_combined
Characteristic
Depression Model
Behavioral Problems Model
OR1
95% CI1
p-value
OR1
95% CI1
p-value
ACEs
0 ACEs
—
—
—
—
1 ACE
1.21
0.99, 1.48
0.067
1.15
1.00, 1.33
0.042
2-3 ACEs
1.38
1.14, 1.66
0.001
1.12
0.97, 1.28
0.11
4+ ACEs
1.70
1.38, 2.09
<0.001
1.37
1.17, 1.59
<0.001
Child Race
Non-Hispanic White
—
—
—
—
Non-Hispanic Black or African American
0.97
0.73, 1.29
0.8
1.04
0.87, 1.26
0.7
Non-Hispanic American Indian or Alaska Native
1.33
0.67, 3.02
0.5
0.93
0.58, 1.58
0.8
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander
1.30
0.80, 2.24
0.3
1.10
0.81, 1.53
0.6
Non-Hispanic Other or Mixed Race
1.01
0.80, 1.30
>0.9
0.92
0.77, 1.10
0.3
Hispanic/Latino of Any Race
0.83
0.69, 1.02
0.067
0.92
0.80, 1.07
0.3
Child Sex
Male
—
—
—
—
Female
1.44
1.26, 1.63
<0.001
1.20
1.08, 1.34
<0.001
Federal Poverty Level
<100%
—
—
—
—
100%-199%
0.88
0.70, 1.11
0.3
0.89
0.75, 1.06
0.2
200%-299%
0.83
0.65, 1.05
0.13
0.81
0.68, 0.97
0.024
300%-399%
0.78
0.60, 1.01
0.059
0.66
0.55, 0.80
<0.001
≥400%
0.67
0.53, 0.85
<0.001
0.69
0.58, 0.83
<0.001
Caregiver Education
Less Than High School
—
—
—
—
High School/Vocational
1.11
0.72, 1.67
0.6
1.00
0.71, 1.37
>0.9
Some College/Associate Degree
1.16
0.75, 1.73
0.5
0.97
0.70, 1.33
0.9
College Degree/Higher
1.05
0.68, 1.57
0.8
0.85
0.61, 1.17
0.3
Child Age
0.96
0.94, 0.98
<0.001
0.86
0.84, 0.87
<0.001
1 OR = Odds Ratio, CI = Confidence Interval
Generate odds ratio plots for Depression and Behavioral Problems stratified by flourishing status. First, create separate plots for the flourishing group and the non-flourishing group, adjusting the x-axis limits, confidence intervals, and dot size. The plots are combined for each group to visualize the odds ratios for the Depression and Behavioral Problems models, and labeled accordingly for each condition within both flourishing and non-flourishing samples.
Expand Code
# Combine the flourishing plots into a list vectorall.flourish.models <-list()all.flourish.models[[1]] <- flourish_model_depall.flourish.models[[2]] <- flourish_model_beh# Plot the combined versionplot_models(all.flourish.models, axis.lim =c(0.2, 5), dot.size =1, show.p =TRUE, ci.lvl =0.95, title ="Odds Ratios by Flourishing Stratified Model", legend.title ="Model", m.labels =c("Depression", "Behavioral Problems"))
Expand Code
# Combine the non-flourishing plots into a list vectorall.non.flourish.models <-list()all.non.flourish.models[[1]] <- non_flourish_model_depall.non.flourish.models[[2]] <- non_flourish_model_beh# Plot the combined versionplot_models(all.flourish.models, axis.lim =c(0.2, 5), dot.size =1, show.p =TRUE, ci.lvl =0.95, title ="Odds Ratios by Non-Flourishing Stratified Model", legend.title ="Model", m.labels =c("Depression", "Behavioral Problems"))
Expand Code
# Create a diagram representing the analysis flowdiagram <-grViz(" digraph flowchart { # Node definitions node [shape = box, style = filled, fillcolor = lightgrey, fontname = Helvetica] X [label = 'Adverse Childhood Experiences', width = 2] Y [label = 'Anxiety, Depression, & Behavioral Problems', width = 2] W [label = 'Child Flourishing', shape = box, style = filled, fillcolor = lightgrey, fontname = Helvetica] Interaction [label = 'ACEs * Child Flourishing Interaction', shape = box, style = filled, fillcolor = lightgrey, fontname = Helvetica] # Define edge connections X -> Y [label = 'a: ##'] W -> Y [label = 'b: ##'] Interaction -> Y [label = 'c: ##'] # Interaction effect X -> Interaction [label = '', style = invis] W -> Interaction [label = '', style = invis] # Graph layout rankdir = LR; # Left to right layout }")# Render the diagramdiagram
Session Information
To enhance reproducibility, details about the working environment used for these analyses can be found below.